
#ifndef MT02
#define MT02

#define EPS 1.0E-5
#define JMAX 30
#define K 5

using namespace std;

static int tcount = 0;

//static double s;

void intPol(double xa[], double ya[], int n, double x, double &y, double &dy){

    int i, m, ns = 0;

    double den, dif, dift, hd, hu, w, *u, *d;

    u = new double[n];
    d = new double[n];

    assert(u != NULL && d != NULL);
    dif = fabs(x - xa[0]);

    for (i = 0; i < n; i++){
        if ((dift = fabs(x - xa[i])) < dif){
            ns = i;
            dif = dift;
        }

        u[i] = ya[i];
        d[i] = ya[i];
    }

    y = ya[ns--];
    //cout << "Y0: " << y << endl;
    for (m = 1; m < n; m++){
        for (i = 0; i < n - m; i++){

            hd = xa[i] - x;
            hu = xa[i + m] - x;

            w = d[i + 1] - u[i];
            den = hd - hu;

            if (den == 0)
                break;

            den = w/den;

            d[i] = hd * den;
            u[i] = hu * den;
        }

        if (2*ns < (n - m - 1))
            dy = d[ns + 1];
        else
            dy = u[ns--];
//        cout << "dY pokytis " << m << " zingsnyje, nuo nario " << ns+2
//             << "(" << y << "): " << dy << endl;
        y += dy;
        //cout << "Y" << m << ": " << y << endl;
    }

    delete []u;
    delete []d;
}



double trapec(double (*func)(double), double a, double b, int n){
	tcount++;
    double x, tnm, sum, del;
    static double s;

    int j, i, it;

    //int it;

    if(n == 1){
        return (s = 0.5 * (b-a)*(func(a) + func(b)));
    }

    for(it = 1, i = 1; i < n - 1; i++)
        it <<= 1;
    tnm = it;
    del = (b - a) / tnm;

    x = (a + 0.5*del);

    for(sum = 0.0, j = 1; j <= it; x += del, j++){
        sum += func(x);
    }

    return (s = 0.5*(s + sum*del));
}


double sumTrapec(double (*func)(double), double a, double b){
    int j;
    double SNew, SOld = -1.0E-30;

    for( j = 1; j < JMAX; j++){
        SNew = trapec(func, a, b, j);

        if(fabs(SOld - SNew) < EPS*fabs(SOld)){
            return SNew;
        }

        if(SNew == 0.0 and SOld == 0.0 and j > 6){
            return 0.0;
            // return SNew;
        }

        SOld = SNew;
    }

    cout << "Per daug iteraciju; rezultato nera" << endl;
    return 0.0;
}

double sumTrapecSyms(double (*func)(double), double a, double b){
    int j;

    double SOld, SNew, simNew, simOld;

    simOld = SOld = 1.0e-30;

    for(j = 1; j < JMAX; j++){
        simNew = trapec(func, a, b, j);

        SNew = (4.0 * simNew - simOld)/3.0;

        if(fabs(SNew - SOld) < EPS * fabs(SOld)){
            return SNew;
        }

        SOld = SNew;

        simOld = simNew;
    }

    cout << "Nieko" << endl;

    return 0.0;
}


double sumTrapecRondberg(double (*func)(double), double a, double b){
    int j;
    double S, SErr;

    double hMas[JMAX], SMas[JMAX];

    hMas[0] = 1.0;

    for( j = 1; j < JMAX; j++){
        SMas[j - 1] = trapec(func, a, b, j);
        if(j >= K){
            intPol(hMas + j - K, SMas + j - K, K, 0.0, S, SErr);
            if(fabs(SErr) <= EPS * fabs(S))
                return S;
        }

        hMas[j] = 0.25 * hMas[j - 1];
    }

    cout << "Daug zingsniu; nieko neturim" << endl;

    return 0.0;
}


#endif

/// Padarykite funkciją integruojančią
// Integruoti funkcijas knygutėje
// y = x^4 * log(x + sqrt(x^2 + 1))
// y = 10 * x, kai 0 < x < 0.5
//   = (20 - 10*x)/3, kai 0.5 <= x < 2
// palyginti našumą


